Collective Langevin Dynamics of Flexible 
Cytoskeletal Fibers 



Francois Nedelecf and Dietrich Foethke 

European Molecular Biology Laboratory 

69117 Heidelberg, Germany 

f To whom correspondence should be addressed 

E-mail: nedelecQembl.de 

Abstract. We develop a numerical method to simulate mechanical objects in a viscous 
medium at a scale where inertia is negligible. Fibers, spheres and other voluminous objects 
are represented with points. Different types of connections are used to link the points 
together and in this way create composite mechanical structures. The motion of such 
structures in a Brownian environment is described by a first-order multivariate Langevin 
equation. We propose a computationally efficient method to integrate the equation, and 
illustrate the applicability of the method to cytoskeletal modeling with several examples. 
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1. Introduction 

The internal arcliitecture of living cells relies largely on microscopic fibers, which form 
the cytoskeleton with their associated proteins. These fibers have remarkable mechanical 
properties. Microtubules and actin filaments for instance have persistence lengths of 
~5 mm and 20 /im, respectively, and can sustain pico-Newtons of force without breaking 
[1]. Yet these fibers can also be broken down quickly, because they are formed by the 
non-covalent assembly of protein monomers. Filament ends can grow or shrink, or even 
alternate between those two states in a remarkable process called dynamic instability [21 13]. 
Structurally, the monomers in microtubules and actin filaments assemble head to tail in a 
regular manner. On the resulting polar lattices, mechano-enzymes called molecular motors 
(for example kinesin on microtubules or myosin on actin-filaments) use chemical energy 
to move directionally or to organize the filaments in space Furthermore, specific 
enzymes control the filaments by regulating nucleation, assembly /disassembly or even by 
severing the filaments. 

The cytoskeleton is involved in multiple cellular processes such as cytokinesis, motility, 
polarization and mitosis. These functions are accomplished by many filaments working 
together. In this way, a set of dynamic or short-lived filaments may form a stable larger 
assembly, as exemplified by the mitotic spindle Many of the enzymes involved in the 
assembly of these structures are part of multi-functional entities El E]. For example, 
motors form oligomers that can actively connect filaments together |1] ; motors may be able 
to disassemble filaments [6]; nucleation can be controlled such that it occurs on existing 
filaments [HI Ej; crosslinkers may be polarity-specific [10] and motors are sometimes linked 
to proteins that track the tips of growing microtubules [3 [TTl [121 US] • Generally speaking, 
modularity allows the cytoskeleton to be reprogrammed, for example at different stages of 
the cell cycle. It allows cells to reuse the same functional elements to achieve different tasks 
and multiplies the number of way in which the organization of fibers can be regulated. This 
modularity is certainly a consequence of the combinatorial exploration operating during 
natural selection [H]. In any case, the cytoskeleton in addition to fibers contains a kit of 
activities which can be combined in many ways. 

Biological systems are hard to understand, and theory is necessary to approach the 
non-intuitive aspects [IS]. It is notable that many models in the cytoskeleton field often 
include the same basic elements (for a recent review on this subject, see [IE]). This 
reflects the inherent modularity of the biological design illustrated briefly in the previous 
paragraph, and also affects the modeling approach. It implies that it is worthwhile to build 
a computer simulation to model a few basic elements, if these elements can be combined 
freely to rapidly model diverse situations. In practice, the elements of the simulation {eg. 
a model of kinesin, or a model of a severing enzyme) can even be implemented, tested 
and benchmarked by different teams of experts for each aspect of the system. Sharing 
computer code in this way can in fact be a practical mean to combine the efforts of the 
community. 
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Writing a cytoskeletal simulation is likely to be a collective task also because it is a 
demanding project, involving multiple aspects: (a) chemical reactions that occur inside 
cells, (b) transport along fibers, for example the motion of molecular motors, (c) assembly 
dynamics of cytoskeletal fibers and (d) motion and deformation of fibers. Fortunately, 
numerous algorithms are available for certain of these aspects, in particular for the reaction- 
diffusion (see [m HE])- Transport along fibers can be modeled with advection equations, 
or with more details of the motion of the motors [19]. The assembly dynamics of fibers has 
been the subject of much research and cannot be reviewed here (see [IS])- The deformation 
of the fibers is a classical mechanical problem (see for example [201I2T] ). However, the scale 
of living cells is associated with many specific features. In particular, Brownian motion 
plays a fundamental role, inertia is negligible [22] and the fibers are dynamic: they can 
lengthen or shorten by self-assembly. As a consequence, the physics of biological fibers 
is fundamentally distinct from other mechanical systems. In brief, public or commercial 
codes are not adapted to simulate the cytoskeleton. 

The purpose of this paper is to describe a method to calculate the mechanics of 
an ensemble of connected fibers and other objects, which is the basis of a cytoskeletal 
simulation such as cytosim. The physics of such system is described by a Langevin 
equation (for an introduction, see [23]) that recreates the Brownian motion of the fibers 
and includes bending elasticity, fiber-fiber interactions and external force-fields. Following 
earlier work [2^ [25] , we use constraints in order to maintain the length of the fibers. This 
is an alternative to methods in which potentials are used to represent the longitudinal 
stiffness of fibers. We extend this approach by introducing an implicit integration scheme. 
Our method was first used to simulate the effects of motor complexes on two radial arrays 
of microtubules (asters) [26], and more recently the assembly of anti-parallel microtubule 
arrays in S. pombe [7] and the positioning of the spindle in the C. elegans embryo [2Z]- A 
major aim of these simulations was to reconstitute the system's operation in silico, from 
established physical principles. This offers two major advantages: i) the assumptions of 
the model are well defined and can always be modified; ii) any property of the system 
can be measured easily. This facilitates further investigations. For example we could 
systematically simplify the model in order to identify a minimal set of working properties 
[7]. In addition, we could identify the parameter range under which the system can 
operate [2Z]. However, for these results to be valid, the systems operation needs to be 
reproduced correctly at the first place! To maximize the chances of success, it is desirable 
to reconstitute the mechanics in a physically sensible and accurate way. One may otherwise 
derive conclusions which do not apply to the real system. 

In this paper, we focus on the mechanical aspects of the fibers, and explore the 
numerical resolution of the associated equations. We first describe objects that in addition 
to fibers are useful for simulating different cellular skeletons. We then present the equation 
of motion and discuss its numerical integration. We examine the numerical stability of the 
resulting method and discuss how it affects the simulation speed. Finally, we discuss how 
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other aspects of the cytoskeleton can be added to extend the mechanical calculation. 

2. Objects 

More accurate mechanics can be achieved if we introduce two new objects in addition 
to fibers: spherical sets of points [spheres) and non-deformable sets of points (solids). 
These objects are also described with points but have different morphologies (see fig. [T]). 
The mechanical properties are also distinct. While fibers may bend, the solids do not 
deform. The spheres can represent spherical viscous membranes such as vesicles. Any 
number of objects can be combined in various ways to build complex cytoskeletons. For 
example, to simulate interacting microtubule asters fibers were positioned around a 
solid using static links (see fig. |2]A.). The solid represented in this case the organelle (called 
the centrosome) which in the cell generates microtubules in a radial fashion. In vivo as well 
as in the simulation, the resulting structure is radially symmetric, and the fibers have their 
ends mechanically joined together. Two such asters were further connected by another 
solid, to model the positioning of the mitotic spindle in C. elegans (2^. In this case, the 
additional solid represented the pole-to-pole mechanical connection achieved by the mitotic 
spindle. To simulate nuclear positioning in S. pombe, fibers (microtubules) were attached 
to a sphere, and the ensemble was confined in a cylindrical volume (see fig. [2^)- The fibers 
and the sphere represented microtubules and the cell nucleus, which are attached also in the 
real cell. To model the formation of anti-parallel microtubule arrays in S. pombe [7], fibers 
where connected by motors and other crosslinkers (see fig. |2p). Using fibers and solids, 
it is also possible to model the segregation of parM plasmids in E. coli (see fig. [2^3), a 
process which depends on actin-like filaments |28]. The objects can naturally be combined 
in many more ways than illustrated here. This enables diverse cellular mechanics to be 
reproduced, and consequently widens the application scope of the method. This freedom 
is intimately linked to the structure of the master equation that will be examined below, 
and to the way it is integrated numerically. 

3. Constrained Langevin Dynamics 

In the simulation, fibers and other objects are described by points. The coordinates of 
the points are collected in a vector x of size Nd, for a system of points in dimension d. 
Following Langevin (for a simple introduction, see |23]) the equation of motion reads: 

dx = nF{x,t)dt + dB{t) (1) 

F(x, t) of size Nd contains the forces acting on the points at time t. It includes object- 
specific forces such as bending elasticity, and all the links between different objects. dB{t) 
of size Nd summarizes the random molecular collisions leading to Brownian motions; it 
is a stochastic non-differentiable function of time. The matrix fi contains the mobility 
coefficients of the object-points, which will be defined later for each object. 
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In addition, certain distances between points inside the objects (|aj — aj\ = Xij) must 
be conserved during the motion. To satisfy these constraints, we perform a step of the 
dynamics in a subspace tangent to the manifold defined by the constraints, and project 
the result on the manifold. The procedure can be explained simply for a point n constrained 
to move at a distance r from a fixed position uq (see fig. |3]). To calculate the motion of n, 
we first write its dynamics in the plane tangent to the sphere at the current position (this 
is the plane allowed by the constraint |n — riol = r). The restricted dynamics is integrated 
implicitly, and the result projected on the sphere to restore the constraint exactly. This 
approach can be generalized as described next. 

4. Numerical integration 

From an initial configuration, the system is calculated by discrete time steps r (see [22] 
for a general discussion on numerical integration). To calculate ^s^t+r from x^, the equation 
([T| is integrated implicitly. We will discuss the advantages of using an implicit rather 
than an explicit integration in section |9| and concentrate here on the practical issues. For 
an implicit integration, we need to express -F(x, t) linearly as x + Gt, where the square 
matrix At contains the stiffness coefficients associated with the interactions, and the vector 
Gt contains the constant forces. This linearization is obtained by summing over all the 
interactions present at time t (see fig. |6]). In our simulations, many of the interactions were 
modeled as harmonic potentials for simplicity, and are therefore already linear. Non-linear 
interactions simply need to be linearized at this point. In particular, the linearization of 
the constraints leads to an orthogonal projection -P(x), which will be defined later for each 
object. To obtain a finite difference scheme for the interval [t, t + r], P and A are used at 
time t, but x is used at t + r (using Xt+r instead of Xt is the basis of implicit integration): 

xt+^ - Xt = Pt[ Tfi{At xt+^ + Gt) + SBt], 

leading to a system of linear equations: 

[ / - rPtfiAt ] i^t+r - = Pt [ TfxiAt Xt + Gt) + 5Bt], (2) 

where At = A{t), Gt = G{t), Pt = P{xt). The "simulated Brownian" 5Bt = ll^^ dB 
is a vector {Pi9t^i}i^[i^Nd\, where 9t^i ~ N(0,1) are Nd independent normally distributed 
numbers (derived from uniformly distributed pseudo-random numbers [29]). The factors 
A ~ '^^^'^ represent the magnitude of the Brownian motion during a lapse of time r. We 
will see later how they are obtained by calibrating the diffusive motion for the objects. 
The equation can be solved to obtain Xf+^, since both the right-hand side and the matrix 
[/ — rPt/iAf] are known. It would be inefficient to invert the matrix, because the system is 
sparse (it only has few non-zero coefficients). This is true of matrix At, as long as objects 
are only connected to few others. This is also true of Pt which is block-diagonal: it has 
one block for each object on the diagonal, but the rest of the coefficients are null. This is 
because the constraints never involve points from different objects, and the projection can 
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thus be done independently for each object. In this situation, it is advantageous to solve 
the linear system using an iterative method j29] . Different iterative solvers are adapted to 
different matrices. Because PtA^ is non-symmetric, we have used the biconjugate gradient 
stabilized (http://www.netlib.org). This method iteratively converges toward the solution 
of the linear system, and can be stopped when the difference with the exact solution is 
below a certain threshold. We set this threshold to ijj min{f5i) , with ijj = 1/10. In this way, 
the numerical error on x remains below 10% of the Brownian motion, and the approximate 
solution of ^ is practically indistinguishable from the real one. In practice, it is wise to 
systematically vary ip and r for each application to check the convergence of the method. 
It is easy to verify, for example, that more stringent values of ip produce the same results. 

Finally, since equation ^ is obtained by linearization, an additional correction is 
necessary to re-establish the constraints. The result of equation ^ is projected back on 
the manifold associated with the constraints [SB]. This introduces corrections which are 
second-order in r. In the following sections, we will call this procedure 'reshaping' the 
objects. We now survey how fibers, spheres and solids are represented in space, their 
mobility coefficients, projection operators and 'reshaping' procedure. The interactions 
between objects (which contribute to At and Gt) will be described subsequently. 



5. Linear set of points (fiber) 

Fibers are modeled as infinitely thin linear objects behaving like elastic, non extensible 
rods [26]. Each fiber is represented by p + 1 equidistant model-points rrii, for i e [0,p], 
separated by a distance L/p. A fiber is polar: mo is the minus-end and nip the plus-end. 
The number of segments p is adjusted as a function of the total length L of the fiber. 
Points are added or removed, in order to always minimize \p — L/p\, for each fiber as it 
grows or shrinks (see fig. |4]). The desired segment length p is a parameter affecting the 
precision of the simulation. To set p, one may run a representative case with various values 
(for microtubules, p < 0.5 pm is usually appropriate). 

It is often necessary to interpolate between the model-points, when for example 
calculating the position x of a molecule attached to the fiber. If and ruk+i are the 
model-points on each side of x, we use x = {l — a)mk + amk+i- The interpolation coefficient 
a G [0, 1] is calculated from the known relative positions of the three points along the fiber: 
= \fnkx\/\mkmk+i\. The model-points are themselves updated using this interpolation 
procedure at every time-step if the length of the fiber has changed (see fig. 111). 



5.1. Bending elasticity 

Fibers can bend under external forces and resist these forces elastically. The standard 
formula for bending elasticity [20] can be applied to strings of points. For any set of three 
consecutive points m^, k & {i — 1; i; i + 1}, we approximate it linearly as a triplet of 
forces {—F; 2F; —F}. Each triplet corresponds to the torque generated between two 
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consecutive segments (see fig. |5]). Furthermore, we have F = a(mj_i — 2mj + mj+i), with 
a = n{p/L)^, where k is the bending modulus of the fiber, and L/n the length of each 
segment. The result was verified by comparing the buckling threshold in the simulation 
with Euler's formula vt^k/L^. The procedure is appropriate if p is such that the angles 
between consecutive segments remain small during the simulation (not shown). Physically, 
the forces are isotropic, i.e. they can be written as a reduced matrix of size p x p (and not 
pcixpc?), obtained by adding several times the 3 X 3 matrix = —(1, —2, 1)®(1, —2, 1) (0 
is the tensor product). The final result is simple because points are distributed regularly 
over the length of the fiber (see fig. [s]). 



5.2. Mobility 

The motion of an object at low Reynolds number is characterized by a mobility. This is 
defined by factors which link speed and force (speed = mobility x force). These factors 
depend on the size and shape of the object, and on the viscosity r] of the surrounding fiuid. 
For instance a straight cylinder has two mobility factors, because it is twofold easier to move 
in the longitudinal direction than in a transverse direction. This anisotropy could not be 
implemented simply, because fibers in the simulation may bend and adopt arbitrary shapes. 
An exact calculation would require finding the hydrodynamic interactions between all the 
points in the system. This can be done in the future, but for simplicity, we have so far used 
the averaged mobility of a straight rod of length L and diameter 6: p = \og{Lh/S)/3nriL 
[50] . The logarithmic term is an effective hydrodynamic correction on the scale L^, which 
is either the length of the fiber, or a hydrodynamic cut-off, whatever is smallest. We derive 
a single mobility factors for the p + 1 points representing a fiber: Pp = {p + 1) /i. 



5.3. Projector associated with the constraints 

In this section, we calculate the projection P derived from the constraint that the length of 
the fiber should remain constant during the resolution of equation ([T]). For each fiber, the 
coordinates of the p + 1 model-points nik are stored in a vector of dimension [p + l)d (for 
d = 3, {xo,Xi,X2} correspond to tuq, and {a;3,X4,a;5} to mi, etc). The motions of these 
points are determined by external forces f = {fk}, and additionally by internal forces 
f = {fk}- The speeds resulting from f + f should be compatible with the constraints 
Ck = {fTLk+i — mkY — {L/pY = for k G [0,]9[. To calculate f from f, we first define the 
p X d{p + 1) Jacobian matrix Jjj = dCi/dxj. In 3D, it reads: 



J = 2 



f X0 — X3 X1—X4 X2 — X5 X^ — Xq X4 — X1 X^ — X2 ■■■\ 

X3—XQ X4 — X7 X5 — XS XQ—X3 Xy — X4 Xs—X^ ■■■ 

'■■ '■■) 



Because the mobility coefficients are the same for all the points (/ip, see sec. 5.2), the 
speed of the points is v = /ip(f + f). This motion maintains the constraints if Jv = 0. 
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Therefore f must be such that J(f + f) = 0. Furthermore, internal forces should not 
contribute to global motion or rotation of the object. This imposes that their work should 
be null for any motion compatible with the constraints: f -u = for any u such that J u = 0. 
This implies that f = J*A, where A is a vector of size p (the Lagrange multipliers). We 
derive J(f + J*A) = 0, and since JJ* of size p x p is non-singular, A = —{JJ^)~^Jf, 
and finally f = —J^{JJ'^)~^Ji. This shows that the total force can be obtained linearly 
as f + f = Pf, with P = I — J^{JJ*)~^J. From this result, it is clear that P is an 
orthogonal projection (P is symmetric and idempotent PP = P). Notice that J J* is 
banded symmetric, and therefore easy to invert, which means that P can be computed 
fast. P (which depends solely on x) is one block of the operator Pt used in equation (j2]). 

Fibers are 'reshaped' to restore the constraints exactly after the model-points have 
been moved. This is done sequentially for k G [0,p], by moving the points ■mQ...mk in the 
direction of mfc+i— mfc and rUk+i-.-mp in the opposite direction, to restore \mk+i—mk\ = L/p 
while conserving the center of gravity of the fiber. 



5.4- Brownian motion 

To simulate Brownian motion, a term 5Bt is attributed to each fiber coordinate Xt (equation 
[2]). This term is most simply calibrated by considering diffusion in the absence of bending 
or external forces (A = 0, G = 0). If we first assume Pj = / in equation (|2]), we get 
Xt+h — Xt = SBt. To produce a pure diffusion with a coefficient D, one needs: 

This holds true if 6Bt is normally distributed, of mean zero and variance 2Dt. We can 
use 6Bt = (36, where 6 ~ N(0,1) is a random number generated for each time step, and 
(3 = V2D^, as mentioned in section [ij From Einstein's relation, we set D = ^pksT, where 
fip is the mobility, ks the Boltzmann constant, and T the absolute temperature. For a fiber 
with p+1 points, we use {p+l)d random numbers, independent and all normally distributed 
of variance Projecting these numbers with P produces the appropriate diffusion for 
the fiber, as well as thermally-driven deformations. For example, the translation x of the 
center of gravity depends on the sum of all the terms in 6B corresponding to the fiber, 
leading to a diffusion D = nksT (with /i and not Hp). 



6. Spherical set of points (sphere) 

To simulate the nucleus of S. pombe and attach microtubules on its surface (see fig. [2|3), 
we implemented a 'spherical set of points' of radius r. Such object is composed of a point 
no in the center, and q additional points Ui on the periphery. If we define = Uk — uq, 
the constraints are \rk\ = r. A sphere moves as a rigid body, and the peripheral points 
behave as if they were embedded in a viscous surface (see fig. fl]). If is the force applied 
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at point k, the motion of the set reads: 

duo = ii^F dt + dB'^ 

/ X / N (3) 

drk = [fi^M dt + dB^) xvk + Pk [fi^ fk dt + dB^) 

where F = J2i=o fi is the total force on the sphere, M = I]?=i '^j ^ /« is the total torque 
calculated from the center, and where 

n 

is the projection on the plane tangent to the sphere in r^. dB^, dB^ and dB^ are the 
Brownian terms. Note that these equations would not describe a set of peripheral points 
articulated around a central node. For example, the motion of the center no depends on 
the sum of all the forces applied to the object, and not only on the force applied in uq. This 
in fact corresponds to a sphere with points on its surface. To keep track of the orientation 
of the sphere, we also included three reference points hk on the surface, which form with 
riQ a reference frame associated to the sphere. The motion of these reference points is 
entirely determined by the total torque on the sphere: df^ = {ji^M dt + dB^^ x f^, where 
as before fk = fik — nQ. When the object needs to be 'reshaped', the peripheral points are 
simply projected on the surface (rio is not moved). 



6.1. Mobility and Brownian Motion 

The equations involve three mobility factors: the translation and rotational mobility of 
the sphere fi^ and fi^, and the mobility of the points in the surface fx^ . Stokes' law 
can be used to set fi^ and fi^, if the sphere is surrounded by a large volume of fluid. 
The mobility coefficients for the points in the surface can also be calculated [21]. As 
described above, points undergo three different types of motion, and a random number SBt 
in equation ^ is associated with each of these motions. The parameters are calculated by 
considering diffusion in the absence of other forces {A = and G = 0). For the translational 
diffusion of the sphere, the result from equation ^ is obtained as previously for the fiber: 
P'^ = ^2fi'^TkBT. Rotational diffusion is calibrated using equation (3). If is fixed on 
the surface, we get rt+r — = SB^ x r^. This should be a rotational diffusion of a point 
on a sphere: 

{rt+r - n) = 0, {{rt+r - n)') = 4 kBTfiVr. 

Since \rt\ = r, we can use for 6Bl^ a random vector with d independent components of 
mean zero and variance 2Tfi^kBT/r'^. A peripheral point rj also diffuses on the surface, 
which in equation ^ is described by rt+r — = Pk SBk,f The projection pt of should 
diffuse in 2D: 

{pt+r - Pt) = 0, {{Pt+T - Ptf) = 4 ksTfi'^T. 

Since Pk is the identity in the tangent plane, we used for SB^^ a vector with d independent 
components of mean zero, and variance 2 r^^ksT. 
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7. Non-deformable set of points (solid) 

We also implemented non-deformable objects called solids (see fig. [T]) in which the points 
move together in such a way that the shape and size of the set is conserved. The number 
of points p in a solid, and their positions Si can be chosen arbitrarily, and each point 
is associated with a radius > 0. The mobility of the solid is derived from Stokes's 
result for the spheres of center Si and radius a^, neglecting for simplicity the hydrodynamic 
interactions between the spheres. It is possible to include points with Oj = provided 
that > ^- In our previous work, we have actually used solids where only one was 
non-zero. These solids moved like isolated spheres, and the points where positions to 
which forces could be applied. 

7.1. Mobility and Constrained Motion 

Because the set of points should not deform, its elementary motion during a time-step can 
be written as {sf^'^ — s\)/t = v + u x s*, where v and uj are instantaneous translation and 
rotation speeds. The spheres of radius Oj in a medium with viscosity rj have a translational 
drag coefficient = Gnrjai, and a rotational drag coefficient = Snija^ [30]. The forces 
and torques resulting from the friction of the fluid on the sphere thus read: 

/, = ^.{v + uxsi), Mi = i^u, (4) 

and should match the externally applied forces fi. 

Y.~fi = Y.fi^ X /i + Mi = X /i (5) 

i i i i 

This set of four equations can be solved algebraically in both 2D and 3D, to express v 
and u; as a function of the external forces /j. The result always flts in the format of 
equation ([T|. It is actually not necessary to calculate the matrix P to run a simulation. 
It is more efficient to calculate v and u when the product P^f is needed. To 'reshape' a 
solid, one may restore a reference conflguration in the current position and orientation. For 
this, the best translation and rotation which brings the reference points onto the current 
points is calculated [32] . The current points are then replaced by the transformed reference 
conflguration. The Brownian components are calibrated as described before. 

8. Interactions between objects 

The three objects deflned previously can be linked together using elementary interactions. 
By adding the contributions of all these interactions in the system, we obtain the linearized 
force F(x, t) = A^x + Gt, which enters equation (|2]). In practice, each elementary 
interaction leads to a small matrix, which needs to be added to the matrix At and vector 
Gt, at the right rows and columns to correspond to the appropriate points (see example 
on flgure [6]). It is necessary to repeat the procedure at every time step, because the 
position of the interactions may change with respect to the model-points. We deflne four 
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interactions in the case where they connect model-points of the objects. We later explain 
the procedure to connect intermediate positions between the model-points. This approach 
can be generalized to more complicated interactions if necessary. For example, it is possible 
to implement a ring able to slide along a fiber with viscous resistance [33] . 

8.1. Connecting an object to a fixed position. 

The simplest way to immobilize an object is to attach a point a within the object to a 
fixed position g. If the stiffness of the link is k, the resulting force is fa = k {g — a). In 
practice, this means adding —k at one diagonal position in matrix At, and kg to the vector 
Gt (see fig. |6|. Such interactions are used to model gliding assays (see fig. [s]) in which 
motors immobilized on a surface propel fibers in solution. Each attached molecular motor 
leads to an elementary interaction where g corresponds to the place of immobilization, and 
a corresponds to the position on the fiber at which the motor domain is attached. 

8.2. Connecting two objects. 

Points from two different objects can be connected by a link of stiffness k. The forces 
between the points are fa = —ft = k{b — a). These elementary interactions are effective 
to model oligomeric motors [22] and more generally any entity able to connect two fibers 
together (see fig. [2p). In the case of an oligomeric motor, a and b are the positions to which 
the two motor domains are attached on the fibers. 

8.3. Confinement in a convex shape. 

To confine the objects inside a convex shape, we use a harmonic potential that is flat inside 
the allowed region, and rises quadratically away from its edge. Hence, a point a outside the 
cell volume is subject to a force /(a) = k{p{a) — a), where p{a) is the closest point to a on 
the edge of the allowed volume. Because p is also the orthogonal projection of a, the force 
corresponds to a friction-less edge. We linearized f as x ^ k [ca ■ {j){o) — x) ) Ca, where Ca 
is a unit vector in the direction of p{a) — a. This linearization corresponds to the tangent 
plane in p{a), and usually gives a good approximation of /(a) as long as the curvature 
is small. To confine a fiber, it is sufficient to follow the procedure for its model-points, if 
the volume is convex, which is the case for example of the cylindrical yeast S. pombe (see 
fig. [2^). To confine the nucleus of radius r in the same volume, we used a cell volume 
reduced by r. In this way only the center of the sphere needs to be tested. 

8.4. Connecting two objects at a given distance. 

A Hookean spring of stiffness k with a non-zero resting length r between two points a and 
b corresponds to: 
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with 6 = a — b. This force should be hnearized for \S\ ~ r, leading for a to a term krS/\S\ 
in Gt and a contribution in At which is: 



k 



6^6 



if 1^1 < r and 



1^1 



[I 



S®6 
52 J 



otherwise, 



(6) 



and the opposite contributions for b. This interaction can be useful to introduce a repulsion 
between the points. It can for example represent the physical interaction between the 
nuclear membrane and the microtubules in S. pombe (see fig. ^ 



8.5. Interpolation of forces 

We have discussed connections which were attached to model-points. However, in the 
case of a fiber, a molecule may bind at any position x, which is likely to be between 
two model-points nik and rrik+i. When this happens, a is interpolated from the flanking 
model-points using a coefficient a = |mfca;|/|mfcmfc_|_i| in [0, 1]. In the same way, a force / 
applied in x can be distributed to the model-points as = (1 — a)f and fk+i = af. Since 
this procedure preserves any linearity in the relationship between force and coordinates, 
the different matrix elements mentioned previously can be used with interpolated points, 
provided they are multiplied left and right by an appropriate weight matrix. We can 
illustrate the procedure for the simplest connection fa = —fb = k{b — a)of stiffness k 
between two points a and b (section [8.2 ), which reads: 





■(1 

When a and b are model-points, this 2x2 matrix is a reduction of A, corresponding 
to the X, y or z- subspaces. This is sufficient in this case because a Hookean spring of 
null resting length is isotropic, that is to say it does not mix x, y and z coordinates, and 
applies similarly to each subspace. This is not the case for all interactions discussed in 
this section, and it is often necessary to calculate a full matrix. Moreover, when a and 
b are intermediate positions between the model points, we have two indices k, I and two 
interpolation coefficients a,P such that a = (1 — a) rrik + amk+i and b = (1— /3) rni+P m^+i. 
If we define a = 1 — a and /3 = 1 — /3, and 



w 




we get: 



w 



( fk \ 

fk+l 
fl 

\ fl+l J 

The resulting 4x4 matrix is w 




-kw^ 



w 



\ 



k) w^, with 



nil 
V m+i 
,—(3,— (3). We derive that 



a matrix made by adding multiple such interactions is symmetric negative-semidefinite 
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{x'^Ax < 0, for any x). The fact that this is true for any configuration of the connections 
guarantees the numerical stabihty of the method, as explained next. 

9. Numerical Stability and Performance 

We have described all the components of equation ^ which describes the collective 
mechanics of cellular fibers and other objects. The necessary steps of the calculation are 
summarized in figure [7} It is useful at this stage to examine the method mathematically. 
This is usually done by looking at two properties: precision and numerical stability 
The precision is a measure of how the typical error behaves when the time-step r becomes 
small. The numerical stability is a measure of how large r can be, before the calculation 
fails. Numerical precision is important for deterministic equations, for example to predict 
the trajectories of celestial bodies. However, this is not so critical at the cellular scale. In 
fact, to simulate the Brownian motion present in the cell, a random term 6B ~ y/r was 
included in equation The presence of this 'noise' indicates that the physics itself limits 
the precision at which the position of an object can be predicted. This fact undermines 
the usefulness of high precision schemes. The implicit method that we have described is 
of order one: the step's error scales like O(r^), which is better than the physical 'noise' in 
^/T. We found that it was not practically useful to use higher order numerical schemes. 

In contrast, the numerical stability of the method is most important. Indeed, explicit 
schemes usually converge only if the time-step is small. In general, a condition like Tfik < 1 
must be fulfilled, where /j, is the mobility of a point in the system, and k the stiffness of 
the interaction potential. For example, we looked at a test-case in which a microtubule is 
pushed by immobilized motors (see [25] and Fig. |8]). It can be simulated explicitly only if 
r < 1 /is, but the implicit method can use larger time-steps. To achieve this stability, we 
treated the repulsive and attractive interactions in the system differently. Compressive 
forces in the fibers (which are repulsive in nature) were replaced by constraints. All 
the other forces were attractive. This ensured that At would be negative-semidefinite 



(this result was proven in section 8.2 for Hookean interactions of null resting length). 
Mathematically, because Pt is an orthogonal projection, we can show that the eigenvalues 
of I — TuPtAt are always greater than 1, for any value of r. This implies that our integration 
scheme is unconditionally stable. For the other elementary interactions, some instabilities 
may appear, but only for very high values of the time step (not shown). 

Beyond stability, other considerations naturally limit the choice of r. In particular 
the iterative solver might not converge when r is large. The optimal time-step generally 
depends on the problem studied, and it is best to perform systematic trials to find it. 
For the test-case (see fig. |8|, the results are consistent for r < 20 ms. This means that 
a value of 5 or 10 ms would be appropriate. The computational requirements depend 
on the total number of steps (total time/time-step), but also on the cost of individual 
steps. An implicit step of integration is always more costly than an exphcit step, because 
a linear system must be solved. However, the use of sparse matrix techniques reduces 
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the additional work. In practice the considerable reduction in the number of steps makes 
implicit simulations faster (in the test-case, this gain is 10^, using r = 10 ms instead of 
1 fis). Increasing the execution speed is essential if many simulations need to be performed. 
Implicit methods require increased numerical labor, of which we have illustrated the main 
difficulties. Using the method described here, we can simulate the examples shown in 
figure 2 B, C & D much faster than real time using one processor (www.cytosim.org). 

10. Other Elements of a Cytoskeletal Simulation 

In addition to mechanics, a cytoskeletal simulation such as cytosim must include additional 
aspects such as the motion of molecular motors, their binding/unbinding dynamics, as well 
as the transitions between growth and shrinkage of dynamic fibers. These processes can be 
modeled most simply by executing small sub-routines after the Brownian mechanics has 
been calculated, because they correspond to independent operations (see fig. [T]). However, 
two particularly important aspects of cytoskeletal physics need to be mentioned. Firstly, 
only in very particular cases can we approximate the system as a well-mixed reactor. 
At least some of the molecules should be spatially resolved. Secondly, the mechanics 
commonly affects the chemistry. For instance the rates of certain key reactions are force- 
dependent. This is the case for the unbinding rates of molecular motors and for their 
stepping rate (see below). Because these elements are essential for modeling the system 
accurately, it will rarely be possible to apply algorithms developed for purely chemical 
systems {eg. the Gillespie algorithms [3l] or even spatially resolved methods [35]) without 
extensive modifications. We can however use simple and robust simulation strategies, as 
illustrated below in the case of molecular motors. 

10.1. Modeling Molecular Motors 

In cytosim, a motor is characterized by a position, when it is not attached, and by a 
pointer to a fiber and a curvilinear abscissa, when it is attached (see fig. |9|. The abscissa 
is the distance, measured along the fiber, between a reference and the attachment position. 
It is necessary to use a reference which is fixed with respect to the physical lattice, because 
the model-points of a fiber are themselves updated as the fiber grows (see fig. |4]). This 
description neatly separates the details of how the mechanics is implemented from the 
routines simulating the motors per se. This means that the interface with the rest of the 
program can be very simple, with only two procedures: step(f) and attach(m). 

10.1.1. Active Motion. The first procedure step(f) simulates the possible actions of a 
bound motor. The argument / is the load of the motor calculated during the collective 
mechanics. The procedure should decide to detach the motor, or to update the abscissa 
a according to a microscopic model for the interval r. For a well characterized motor 
like kinesin, a classical model is based on the measured characteristics of the motion: the 
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abscissa is increased by Sa = rv^naxi^ ~ ///stall)- In addition, a force-dependent unbinding 
rate pos = Poexp(|/|//o) is used to model the dissociation from the fiber, fmax, Po, fo 
and /stall are characteristics of the motor that have been measured for kinesin With 
this model, the fibers are continuous tracks along which motors may be located anywhere. 
Alternatively, we may model the motion of a motor as a succession of discrete stochastic 
steps. In this case, the motor does one of four things: stay immobile, detach, take a step 
toward the minus-end or take a step toward the plus-end. This means that if the motor 
does not detach, the abscissa is either unchanged, or it is increased or decreased by the step 
size (8 nm). The procedure step(f) calculates the probabilities of these events as a function 
of the force / for the interval r, and selects one of them. This model is quite attractive, 
because these probabilities are actually available for kinesin [36]. Most models describing 
the movement of motors [19] can be summarized similarly with a function step(f). 



10.1.2. Attachment to Fibers. The second procedure necessary to model motors, 
attach (m) simply decides if a unbound motor binds or not to a site m. Usually the 
model would specify e, a maximum distance at which a motor may bind from its current 
position (see fig. |9]). In addition, the molecule would bind at the closest site on the fiber 
(the orthogonal projection) with a certain molecular binding rate fcon (-s"^)- To simulate 
attachments, one therefore needs to first find the fiber-segments which are closer than 
e, typically from all the positions x at which molecular motors are located. For each 
point X, the list of candidates should then be shuffled, to ensure a random ordering of the 
segments. The molecular binding rate can finally be tested sequentially for each segment 
in the list, for example by comparing rkon with a random number 6 in [0,1]. The first 
successful trial is followed by attachment. If done naively, the first step of the operation 
may require calculating the distance of all points to all fiber-segments, and thus a great deal 
of computation for many motors. To avoid this bottleneck in cytosim, a divide and conquer 
algorithm was developed (see fig. 10). Its goal is to limit the number of segments that need 
to be tested to find those which are close to x. The geometrical distance between x and 
these segments is calculated using the vector cross-product to exactly determine which 
ones are closer than e. Reducing the number of tested segments is sufficient to accelerate 
the simulation. 



11. Conclusion 



The method described here is efficient to simulate sparsely connected networks of filaments. 
It applies to many in vivo situations, because the connections between fibers are usually 
mediated by proteins that are small compared to the fibers, and consequently the fibers 
are only locally connected. We have modeled fibers as oriented lines, which is sufficient to 
calculate the extent of bending. It may be necessary in the future to include more details 
such as writhe, since cytoskeletal fibers also have a torsional rigidity. The method can be 
extended in several other ways. One could for instance easily model discrete binding sites on 
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the fibers. This may be important if the fibers are highly covered and molecules compete or 
interact while bound to the lattice. It is also possible to extend the overdamped mechanics 
by adding hydrodynamic effects. It will be very exciting to integrate fiber mechanics 
with membrane dynamics, since membranes and cytoskeleton contribute synergistically to 
cellular architecture, but this might take some time. Cellular chemistry, reaction-diffusion 
of the components in the cell, gene expression networks, can be added more simply. This 
can be done by interfacing our software with other tools (eg. the Virtual Cell project), which 
already cover some of these aspects of physiology. We did not discuss here implementation 
issues, but the scale of the task should remind us of their importance. Software modularity 
is essential to divide the development effort in separate projects of manageable size. Sub- 
models or algorithms should be developed and tested separately, in such a way that they 
can be added or removed from the integrative software easily. Dividing the work among 
different groups is the best way to produce the high-quality cellular simulations that biology 
needs. 
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Fiber Sphere Solid 



Figure 1. Elementary Objects. All objects in the simulation are described by points. 
The points can move in the viscous medium, but the relative distances between certain 
points arc conserved (lines). Left: A fiber is modeled as an equidistributcd string of 
points. Center: A sphere is composed of a central point and peripheral points, located 
a distance r from the center. The peripheral points can move on the surface, as if they 
were in a viscous membrane. Right: A solid is a set of points that behaves like a solid 
body. Its shape and size are constants. 
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Figure 2. Some problems studied with cytosim. In all the images, fibers are 
indicated in white, along with their model-points. (A) An aster is constructed by 
assembling fibers radially around a solid [SSJEl]- Right, top: interactions of microtubules 
with the cell cortex. Right, bottom: the solid is made of a central point (blue) surrounded 
by two concentric layers of peripheral points (green and red). Only the central point is 
associated with a viscous drag (a^ > 0). The other points are used to attach fibers: the 
minus-end to one green point, and a distal position on the fiber to one red point. Using 
a similar simulation with two asters linked by a solid spindle, we proposed an original 
model describing the 3D motions of the spindle in the first cell division of the C. elegans 
embryo [27 . (B) Microtubules in interphase fission yeast and the nucleus, represented by 
a sphere (blue/green). This can be used to study the role of mechanics in regulating the 
dynamics and organization of microtubules. (C) Self-assembly of interphase microtubules 
arrays in fission yeast. The simulation contains no steric interaction between the fibers, 
and they overlap freely. In the display, however, the fibers are shifted in order to visualize 
the bridging complexes (bottom and right). Using this simulation, we could identify a 
minimal 'recipe' to make stable bundles from dynamic microtubules. This recipe describes 
how cross-linking, nucleating and motor activities can be associated to obtain the result 
observed in vivo. (D) Self-segregation of plasmids in prokaryotes. Actin-like filaments are 
simulated, together with two solids, representing the plasmids [2S]. The efficiency of the 
segregation is recapitulated in the simulation, and can therefore be analyzed. 
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Figure 3. Dynamics with Constraints. The principle of the algorithm is illustrated 
here for a point rii constrained to stay at a fixed distance from ng. The point is first 
moved on the tangent to the circle (this is the plane associated with the constraint) using 
an implicit integration scheme. The constraint is then re-established exactly by projecting 
on the sphere. We call this last operation 'reshaping an object'. 




Figure 4. Dynamic Fibers. Top: The model-points of a fiber are updated when 
the tips grow, but they are always equally distributed over the fiber. Points are added 
or removed as necessary to ensure an optimal coverage (see sec. |5|. Bottom: An 
intermediate position x along the fiber is interpolated from the model-points located on 
each side: x — {1 — a)mi^ + am^+i (see sec. jsj. 
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Figure 5. Matrix elements associated with bending elasticity. The stiffness 
matrix At contains the bending elasticity of fibers. The contributions are obtained by 
adding a 3 x 3 elementary matrix for each consecutive triplets of points (see section 5.11. 
The sum of all rows and columns is zero, since the matrix should only generate an internal 
torque. The forces associated with the first triplet (points xi, X2 and X3) are depicted. 
The resulting matrix for 5 points is also shown, and the generalization is straightforward. 
For any fiber, the result is a symmetric banded matrix multiplied by a scalar a that 
depends on the bending elasticity modulus and on the distance between the points. 
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Figure 6. Stiffness matrix and force vector. The stiffness matrix At and the force 
vector Gt in equation [2] are set by considering all the interactions present at time t. For 
each interaction, the appropriate formula (sec. [s]) is first expanded algebraically. The 
factors associated with the coordinates of the points are added to A, and the coefficients 
which are independent of the coordinates are added to G. At the end of the procedure, 
one obtains a (sparse) symmetric matrix A and a vector G that provide the forces on 



the points F = Ax + G. Here we illustrate how a connection of stiffness fci (sec. 8.2 1 



contribute to factors fci and —ki at the rows and columns of A corresponding to the points 



connected. For a connection to a fixed position g (sec. 8.11, a stiffness coefficient — fc2 
is added in A, while is added in G. In this example, the connections are attached 
exactly to points of the system, but this is not always the case. Section [8] explains the 
general proceduce. In addition, the matrix represented here corresponds to a ID system. 



It needs to be duplicated for a 2D simulation, and triplicated in 3D (sec. 8.51 
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Pool coordinates of object-points. 
Calculate projection P for each object. 



Loop over all interactions to set matrix A and vector G. 



Set Brownian components from random numbers, 
record Brownian magnitude in (3. 



Calculate right-hand side of equation (2), 
Solve system of linear equations using iterative method, 

with a precision exceeding (3^, with ^J=0.1. 



Project solution to 'reshape' the objects. 
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For mobile attachments such as molecular motors: 

Calculate tensions in the interaction link. 
Use this information to move attachment positions, 
according to the characteristics of the motors. 



Attachment trials for unbound motors. 
Detachment trials for bound motors, 
(detachment rates are usually force-dependent) 



Calculate forces on fiber tips. 
Elongate fibers according to their force-growth curve. 
Recalculate the model-points of fibers by interpolation. 



Delete/nucleate filaments, add/remove objects. 



Figure 7. Synopsis of a simulation time-step. Sub-steps necessary to simulate a 
system of molecular motors and dynamic fibers. The collective mechanics corresponds to 
the algorithm described in the article. As a byproduct of calculating the mechanics, one 

gets the tensions in the fibers and the forces connecting the fibers. With this information, 
simulation sub-steps can be performed for the objects independently. Events such as the 
binding and the unbinding of motors and the nucleation of new filaments will most likely 
be modeled stochastically. Depending on the level of details required, less-discrctc events 
may be simulated in a deterministic manner. For example, the active motion of molecular 
motors and the assembly dynamics of cytoskeletal fibers can be simulated as non-random 
processes characterized by a force-velocity curve. 
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Figure 8. Numerical stability of the integration scheme. Top: A gliding assay 
where a filament is attached at its end (time-intervals of 5s). The motors pushing the 
fiber lead to the formation of a rotating spiral, as observed experimentally |25j . The 
rotation speed and maximum radius of the spiral can be calculated from the parameters 
of the system: 16000 motors cover an area of 2 x 2/xm, and have the characteristics 
of kinesin: stall force /max — 5pN, unloaded speed OA^m/s, binding rate 10 s"^, 
unbinding rate PofE = 0.5 s^^ exp(force/2.5piV), maximum binding distance lOnm and 
stiffness 200pN/ fj,m. The microtubule of length 8 fim has a rigidity of 20 pN nm^ . It is 
constrained at the minus end by a link of stiffness 4000 pN/ fim. The effective viscosity 
is 0.02 pN s iim^^ . Bottom: The configuration is simulated for different values of the 
time-step r, with accurate results for t < 20 ms. The algorithm is numerically stable, 
and even produces a spiral with r 0.5 s. However, the radius is then under-estimated, 
and the rotation speed overestimated. Another critical parameter, the distance p between 
the points on the fiber was also varied. The results shown for p = 0.1, 0.2, 0.4 and 0.5 
pm (different lines) are similar, because all these values are appropriate. The calculations 
were inaccurate however with p = 0.8 pm (data not shown). This is expected considering 
that the radius of the spiral is ~ I. Apia. 
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Figure 9. Molecular Motors. Top: An unbound motor (diamond) is represented by 
a position. Attachment occurs on the closest site on the fiber-segment, provided this site 
is within a distance e (dashed lines). The capture regions of the segments are truncated 
such that they cover cxac;tly the region located at a distance e from a .straight fiber. When 
the fiber is not straight, the gaps and overlaps exactly compensate each other. Bottom: 
A bound motor is represented by a pointer to a fiber, and by a curvilinear abscissa a{t) 
measured from a fixed origin on the fiber. This defines the position of the motor along the 
fiber independently of the mathematical representation of the fiber. The motor sub-model 
needs to decide whether the motor should detach during the interval of time r, or it needs 
to calculate the displacement da during the same interval. For this, it can use the load / 
calculated during the collective mechanics, and other properties associated with the fiber, 
such as the proximity of the ends, or information on the crowdedness of the binding sites 
on the fiber. 
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Figure 10. Divide and Conquer algorithm. To simulate the attachments to fibers, 
we must be able to find all the fiber-segments which are within a distance e from an 
arbitrary position x. We can proceed according to the following two steps method: Divide 
(left): A grid is set in space, each node of the grid being associated with a list of segments. 
The segments are recorded on the grid, at the nodes located at a distance h or less {h will be 
defined later). This operation is performed in 2D using standard rasterizer codes derived 
from computer graphics, which are optimized to scan all points with integer coordinates 
located inside an arbitrary polygon. We rasterize the rectangles built around the segments 
at a distance h. For example, on this diagram, the blue segment is recorded at the blue 
points, and the red segment at the red points. In 3D, the rectangular volume can be 
rasterized following the same principles as in 2D. Conquer (right): After the segments 
have been distributed over the grid, one can quickly find which ones are near x: one needs 
to check only the segments recorded at the grid point g closest to x. One will find all 
segments located at distance h — d/2 or less from x, since \gx\ < d/2, where d is the 
diagonal of the grid. Hence to find all the segments closer than e, one sets h = e + d/2 
during the rasterizing operation. Note: The grid does not need to be square (the unit 
cell can be rectangular) and it can be adjusted for optimal performance. If the grid is too 
fine, it will use a lot of memory, and rasterizing will be slow. If the grid is coarse (d large), 
the number of candidates returned for a point x will be larger. Experimentation may be 
necessary to optimize the grid, but the procedure provides exact results for any cell size. 
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